!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine model_initial (base)
!==================================================================================================================================|

   use parameters
   use block_define
   use interface_hilbert_curve

   implicit none

   type(block), pointer, intent(inout) :: base

   integer(4) :: i_block, j_block, k_block, i_block_n, j_block_n, k_block_n, i_block_n_p, j_block_n_p, k_block_n_p
   integer(4) :: i_neigh, j_neigh, k_neigh, i_total, j_total, k_total, i_bnd_intra, i_cpu, i_hilbert
   integer(4) :: zerox, zeroy, zeroz
   integer(4) :: average_blocks, blocks_p1_cpu_number, blocks_front, tmp_start, tmp_end

   type(block), pointer :: child, current, next

   real(8) :: pr, b2, v2, c2, s2, cax, cay, caz, cfx, cfy, cfz, vx, vy, vz, x_c, y_c, z_c

!----------------------------------------------------------------------------------------------------------------------------------|
   zerox = nint ((nx_total - 2 * ng) * zero_posx / rangex) + 1 + ng
   base%dx = rangex / dble(nx_total - 2 * ng)

   base%x(zerox) = base%dx * 0.5d0

   do i = zerox + 1, nx_total
      base%x(i) = base%x(i - 1) + base%dx
   enddo
   do i = zerox - 1, 1, -1
      base%x(i) = base%x(i + 1) - base%dx
   enddo

   zeroy = nint ((ny_total - 2 * ng) * zero_posy / rangey) + 1 + ng
   base%dy = rangey / dble(ny_total - 2 * ng)

   base%y(zeroy) = base%dy * 0.5d0

   do j = zeroy + 1, ny_total
      base%y(j) = base%y(j - 1) + base%dy
   enddo
   do j = zeroy - 1, 1, -1
      base%y(j) = base%y(j + 1) - base%dy
   enddo

   zeroz = nint ((nz_total - 2 * ng) * zero_posz / rangez) + 1 + ng
   base%dz = rangez / dble(nz_total - 2 * ng)

   base%z(zeroz) = base%dz * 0.5d0

   do k = zeroz + 1, nz_total
      base%z(k) = base%z(k - 1) + base%dz
   enddo
   do k = zeroz - 1, 1, -1
      base%z(k) = base%z(k + 1) - base%dz
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  initial blocks
!----------------------------------------------------------------------------------------------------------------------------------|
   base%lv = 0
   base%cpu_id = cpu_id
   base%nc_x = ini_block_x
   base%nc_y = ini_block_y
   base%nc_z = ini_block_z
   base%nc = base%nc_x * base%nc_y * base%nc_z
   base%nx = nx_total
   base%ny = ny_total
   base%nz = nz_total
   base%dt = min_dt
   base%t = 0.0d0
   base%var_allocate_flag = .true.
   base%p_pos = (/ 0, 0, 0 /)
   base%id = 0

   if (ini_block_x .le. 1 .or. ini_block_y .le. 1 .or. ini_block_z .le. 1) then
      write (*, *)
      write (*, '(a)') ' Error: ini_block_x or ini_block_y or ini_block_z should be great than 1.'
      write (*, '(a)') ' The program will be terminated.'
      write (*, *)
      call mpi_abort (mpi_comm_world, 1, error)
   endif

   do k_block = 1, ini_block_z
   do j_block = 1, ini_block_y
   do i_block = 1, ini_block_x
      allocate(base%child(i_block, j_block, k_block)%pblock)
      child => base%child(i_block, j_block, k_block)%pblock
      allocate(child%x(nx_block), child%y(ny_block), child%z(nz_block),                                                            &
               child%p_pos(n_dims), child%parent_cell(n_dims, n_cell), child%id(n_dims, n_levels))

      child%lv = 1
      child%p_pos(1) = i_block
      child%p_pos(2) = j_block
      child%p_pos(3) = k_block
      child%id = 0
      child%id(:, child%lv) = child%p_pos
      child%nx = nx_block
      child%ny = ny_block
      child%nz = nz_block
      child%dx = base%dx
      child%dy = base%dy
      child%dz = base%dz
      child%dt = base%dt
      child%t = base%t
      child%parent => base

      do i = 1, child%nx
         i_total = (i_block - 1) * nx_no_bnd + i
         child%x(i) = base%x(i_total)
      enddo

      do j = 1, child%ny
         j_total = (j_block - 1) * ny_no_bnd + j
         child%y(j) = base%y(j_total)
      enddo

      do k = 1, child%nz
         k_total = (k_block - 1) * nz_no_bnd + k
         child%z(k) = base%z(k_total)
      enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!     For the boundary inside the computational box
!----------------------------------------------------------------------------------------------------------------------------------|
      x_c = (minval(child%x) + maxval(child%x)) / 2.0d0
      y_c = (minval(child%y) + maxval(child%y)) / 2.0d0
      z_c = (minval(child%z) + maxval(child%z)) / 2.0d0
      do i_bnd_intra = 1, n_bnd_interior
         if (x_c .ge. bnd_interior_pos(1, i_bnd_intra) .and.                                                                       &
             y_c .ge. bnd_interior_pos(2, i_bnd_intra) .and.                                                                       &
             z_c .ge. bnd_interior_pos(3, i_bnd_intra) .and.                                                                       &
             x_c .le. bnd_interior_pos(4, i_bnd_intra) .and.                                                                       &
             y_c .le. bnd_interior_pos(5, i_bnd_intra) .and.                                                                       &
             z_c .le. bnd_interior_pos(6, i_bnd_intra)) then
            child%lv = -i_bnd_intra
            n_blocks_all = n_blocks_all - 1
            exit
         endif
      enddo

   enddo
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  hilbert space filling curve
!----------------------------------------------------------------------------------------------------------------------------------|
   allocate(unit_cell(n_dims, n_cell))
   unit_cell(:, 1) = (/ 0, 0, 0 /)
   unit_cell(:, 2) = (/ 0, 0, 1 /)
   unit_cell(:, 3) = (/ 1, 0, 1 /)
   unit_cell(:, 4) = (/ 1, 0, 0 /)
   unit_cell(:, 5) = (/ 1, 1, 0 /)
   unit_cell(:, 6) = (/ 1, 1, 1 /)
   unit_cell(:, 7) = (/ 0, 1, 1 /)
   unit_cell(:, 8) = (/ 0, 1, 0 /)

   current => base
   call hilbert_curve(base, current, unit_cell, 3, ini_block_order, (/1, 1, 1/), cpu_id)
   base%last => current
   current%next => base
   base%framework_last => current
   current%framework_next => base

!----------------------------------------------------------------------------------------------------------------------------------|
!  delete the node of internal boundary blocks (lv eq -1) and update value of hilbert_number of every block
!----------------------------------------------------------------------------------------------------------------------------------|
   i_hilbert = 0
   base%hilbert_number = i_hilbert

   current => base
   next => base%next
   i_hilbert = 1
   do while (next%lv .ne. 0)
      if (next%lv .gt. 0) then
         next%last => current
         current%next => next
         next%framework_last => current
         current%framework_next => next
         next%hilbert_number = i_hilbert
         i_hilbert = i_hilbert + 1
         current => next
      endif
      next => next%next
   enddo
   base%last => current
   current%next => base
   base%framework_last => current
   current%framework_next => base

!----------------------------------------------------------------------------------------------------------------------------------|
!  framework gene
!----------------------------------------------------------------------------------------------------------------------------------|
   n_gene = ini_block_x * ini_block_y * ini_block_z

   if (n_gene .le. 0) then
      write(*, *) 'ERROR: Too many levels, the type of integer*4 is overflow.'
      call mpi_abort (mpi_comm_world, 1, error)
   endif

   allocate(i_gene(cpu_number))
   allocate(gene(n_dims, n_levels, n_gene))

!----------------------------------------------------------------------------------------------------------------------------------|
!  allocatememory for exchange data
!----------------------------------------------------------------------------------------------------------------------------------|
   allocate(n_blocks(cpu_number), n_blocks_vip(cpu_number))
   allocate(level_len(n_levels))
   allocate(hilbert_start(cpu_number), hilbert_end(cpu_number))

   allocate(n_blocks_update_send(cpu_number), n_blocks_update_sendrecv(cpu_number, cpu_number))

   allocate(blocks_update_info_send(n_d_l, ratio2 * n_levels, cpu_number))
   allocate(blocks_update_info_recv(n_d_l, ratio2 * n_levels, cpu_number))

   allocate(n_neighs_reflux_sendrecv_x(cpu_number, cpu_number))
   allocate(n_neighs_reflux_sendrecv_y(cpu_number, cpu_number))
   allocate(n_neighs_reflux_sendrecv_z(cpu_number, cpu_number))
   allocate(n_neighs_shareblock_sendrecv_x(cpu_number, cpu_number))
   allocate(n_neighs_shareblock_sendrecv_y(cpu_number, cpu_number))
   allocate(n_neighs_shareblock_sendrecv_z(cpu_number, cpu_number))
   allocate(n_neighs_exchange_sendrecv_x(cpu_number, cpu_number))
   allocate(n_neighs_exchange_sendrecv_y(cpu_number, cpu_number))
   allocate(n_neighs_exchange_sendrecv_z(cpu_number, cpu_number))
!   allocate(n_neighs_exchange_sendrecv_fix_y(cpu_number, cpu_number))
!   allocate(n_neighs_exchange_sendrecv_fix_z(cpu_number, cpu_number))

   reflux_size_x = n_gene
   reflux_size_y = n_gene
   reflux_size_z = n_gene
   allocate(neighs_reflux_info_send_x(n_d_l + 1, reflux_size_x, cpu_number))
   allocate(neighs_reflux_info_send_y(n_d_l + 1, reflux_size_y, cpu_number))
   allocate(neighs_reflux_info_send_z(n_d_l + 1, reflux_size_z, cpu_number))
   allocate(neighs_reflux_info_recv_x(n_d_l + 1, reflux_size_x, cpu_number))
   allocate(neighs_reflux_info_recv_y(n_d_l + 1, reflux_size_y, cpu_number))
   allocate(neighs_reflux_info_recv_z(n_d_l + 1, reflux_size_z, cpu_number))

   shareblock_size_x = n_gene
   shareblock_size_y = n_gene
   shareblock_size_z = n_gene
   allocate(neighs_shareblock_info_send_x(n_d_l + 1, shareblock_size_x, cpu_number))
   allocate(neighs_shareblock_info_send_y(n_d_l + 1, shareblock_size_y, cpu_number))
   allocate(neighs_shareblock_info_send_z(n_d_l + 1, shareblock_size_z, cpu_number))
   allocate(neighs_shareblock_info_recv_x(n_d_l + 1, shareblock_size_x, cpu_number))
   allocate(neighs_shareblock_info_recv_y(n_d_l + 1, shareblock_size_y, cpu_number))
   allocate(neighs_shareblock_info_recv_z(n_d_l + 1, shareblock_size_z, cpu_number))

   exchange_size_x = n_gene
   exchange_size_y = n_gene
   exchange_size_z = n_gene
   allocate(neighs_exchange_info_send_x(n_d_l + 1, exchange_size_x, cpu_number))
   allocate(neighs_exchange_info_send_y(n_d_l + 1, exchange_size_y, cpu_number))
   allocate(neighs_exchange_info_send_z(n_d_l + 1, exchange_size_z, cpu_number))
   allocate(neighs_exchange_info_recv_x(n_d_l + 1, exchange_size_x, cpu_number))
   allocate(neighs_exchange_info_recv_y(n_d_l + 1, exchange_size_y, cpu_number))
   allocate(neighs_exchange_info_recv_z(n_d_l + 1, exchange_size_z, cpu_number))
!   allocate(neighs_exchange_info_send_fix_y(n_d_l + 1, exchange_size_y, cpu_number))
!   allocate(neighs_exchange_info_send_fix_z(n_d_l + 1, exchange_size_z, cpu_number))
!   allocate(neighs_exchange_info_recv_fix_y(n_d_l + 1, exchange_size_y, cpu_number))
!   allocate(neighs_exchange_info_recv_fix_z(n_d_l + 1, exchange_size_z, cpu_number))

   reflux_size_send_x = 1
   reflux_size_recv_x = 1
   reflux_size_send_y = 1
   reflux_size_recv_y = 1
   reflux_size_send_z = 1
   reflux_size_recv_z = 1
   allocate(neighs_reflux_send_x(sys, ny_no_bnd, nz_no_bnd, reflux_size_send_x))
   allocate(neighs_reflux_send_y(sys, nx_no_bnd, nz_no_bnd, reflux_size_send_y))
   allocate(neighs_reflux_send_z(sys, nx_no_bnd, ny_no_bnd, reflux_size_send_z))
   allocate(neighs_reflux_recv_x(sys, ny_no_bnd, nz_no_bnd, reflux_size_recv_x))
   allocate(neighs_reflux_recv_y(sys, nx_no_bnd, nz_no_bnd, reflux_size_recv_y))
   allocate(neighs_reflux_recv_z(sys, nx_no_bnd, ny_no_bnd, reflux_size_recv_z))

   shareblock_size_send_x = 1
   shareblock_size_recv_x = 1
   shareblock_size_send_y = 1
   shareblock_size_recv_y = 1
   shareblock_size_send_z = 1
   shareblock_size_recv_z = 1
   allocate(neighs_shareblock_send_x(sysp1, ng + 1, ny_block, nz_block, shareblock_size_send_x))
   allocate(neighs_shareblock_send_y(sysp1, nx_block, ng + 1, nz_block, shareblock_size_send_y))
   allocate(neighs_shareblock_send_z(sysp1, nx_block, ny_block, ng + 1, shareblock_size_send_z))
   allocate(neighs_shareblock_recv_x(sysp1, ng + 1, ny_block, nz_block, shareblock_size_recv_x))
   allocate(neighs_shareblock_recv_y(sysp1, nx_block, ng + 1, nz_block, shareblock_size_recv_y))
   allocate(neighs_shareblock_recv_z(sysp1, nx_block, ny_block, ng + 1, shareblock_size_recv_z))

   exchange_size_send_x = 1
   exchange_size_recv_x = 1
   exchange_size_send_y = 1
   exchange_size_recv_y = 1
   exchange_size_send_z = 1
   exchange_size_recv_z = 1
   allocate(neighs_exchange_send_x(sysp1, ng, ny_no_bnd, nz_no_bnd, exchange_size_send_x))
   allocate(neighs_exchange_send_y(sysp1, nx_block, ng, nz_no_bnd, exchange_size_send_y))
   allocate(neighs_exchange_send_z(sysp1, nx_block, ny_block, ng, exchange_size_send_z))
   allocate(neighs_exchange_recv_x(sysp1, ng, ny_no_bnd, nz_no_bnd, exchange_size_recv_x))
   allocate(neighs_exchange_recv_y(sysp1, nx_block, ng, nz_no_bnd, exchange_size_recv_y))
   allocate(neighs_exchange_recv_z(sysp1, nx_block, ny_block, ng, exchange_size_recv_z))

!!   allocate(neighs_shareblock_send_x_hc(8, ng + 1, ny_block, shareblock_size_x, cpu_number))
!!   allocate(neighs_shareblock_send_y_hc(8, nx_block, ng + 1, shareblock_size_x, cpu_number))
!!   allocate(neighs_shareblock_recv_x_hc(8, ng + 1, ny_block, shareblock_size_x, cpu_number))
!!   allocate(neighs_shareblock_recv_y_hc(8, nx_block, ng + 1, shareblock_size_x, cpu_number))

!   allocate(neighs_shareblock_send_x_hc(1, ng + 1, ny_block, nz_block, shareblock_size_x, cpu_number))
!   allocate(neighs_shareblock_send_y_hc(1, nx_block, ng + 1, nz_block, shareblock_size_y, cpu_number))
!   allocate(neighs_shareblock_send_z_hc(1, nx_block, ny_block, ng + 1, shareblock_size_z, cpu_number))
!   allocate(neighs_shareblock_recv_x_hc(1, ng + 1, ny_block, nz_block, shareblock_size_x, cpu_number))
!   allocate(neighs_shareblock_recv_y_hc(1, nx_block, ng + 1, nz_block, shareblock_size_y, cpu_number))
!   allocate(neighs_shareblock_recv_z_hc(1, nx_block, ny_block, ng + 1, shareblock_size_z, cpu_number))

!   allocate(neighs_exchange_send_x_hc(1, ng, ny_no_bnd, nz_no_bnd, exchange_size_x, cpu_number))
!   allocate(neighs_exchange_send_y_hc(1, nx_block, ng, nz_no_bnd, exchange_size_y, cpu_number))
!   allocate(neighs_exchange_send_z_hc(1, nx_block, ny_block, ng, exchange_size_z, cpu_number))
!   allocate(neighs_exchange_recv_x_hc(1, ng, ny_no_bnd, nz_no_bnd, exchange_size_x, cpu_number))
!   allocate(neighs_exchange_recv_y_hc(1, nx_block, ng, nz_no_bnd, exchange_size_y, cpu_number))
!   allocate(neighs_exchange_recv_z_hc(1, nx_block, ny_block, ng, exchange_size_z, cpu_number))

!----------------------------------------------------------------------------------------------------------------------------------|
!  update the information of the number of blocks of every cpu
!----------------------------------------------------------------------------------------------------------------------------------|
   average_blocks = n_blocks_all / cpu_number
   blocks_p1_cpu_number = mod(n_blocks_all, cpu_number)
   if (cpu_id .lt. blocks_p1_cpu_number) average_blocks = average_blocks + 1

   call mpi_allgather(average_blocks, 1, mpi_integer4, n_blocks, 1, mpi_integer4, mpi_comm_world, error)

!----------------------------------------------------------------------------------------------------------------------------------|
!  update the value of cpu_id of every block and deallocatethe variables
!----------------------------------------------------------------------------------------------------------------------------------|
   i_block = 1
   i_cpu = 0
   base%cpu_id = cpu_id
   current => base%next
   do while (current%lv .ne. 0)
      if (i_block .gt. sum(n_blocks(1 : i_cpu + 1))) then
         i_cpu = i_cpu + 1
      endif
      current%cpu_id = i_cpu
      i_block = i_block + 1
      current => current%next
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  block boundary
!----------------------------------------------------------------------------------------------------------------------------------|
   do k_block = 1, ini_block_z
   do j_block = 1, ini_block_y
   do i_block = 1, ini_block_x

      child => base%child(i_block, j_block, k_block)%pblock
      do k_neigh = -1, 1
         k_block_n = k_block + k_neigh
         do j_neigh = -1, 1
            j_block_n = j_block + j_neigh
            do i_neigh = -1, 1
               i_block_n = i_block + i_neigh
               child%neigh(i_neigh, j_neigh, k_neigh)%pblock => null ()
               if (i_block_n .ge. 1 .and. i_block_n .le. ini_block_x .and.                                                         &
                   j_block_n .ge. 1 .and. j_block_n .le. ini_block_y .and.                                                         &
                   k_block_n .ge. 1 .and. k_block_n .le. ini_block_z) then
                  child%neigh(i_neigh, j_neigh, k_neigh)%pblock => base%child(i_block_n, j_block_n, k_block_n)%pblock
               else
!----------------------------------------------------------------------------------------------------------------------------------|
!                 periodic boundary
!----------------------------------------------------------------------------------------------------------------------------------|
                  i_block_n_p = i_block_n
                  if (any (bnd_exterior(:, :, 1) .eq. 0) .and. (i_block_n .eq. 0 .or. i_block_n .eq. ini_block_x + 1)) then
                     if (i_block_n .eq. 0) then
                        i_block_n_p = ini_block_x
                     elseif (i_block_n .eq. ini_block_x + 1) then
                        i_block_n_p = 1
                     else
                        write (*, *) 'Periodic boundary error: model_initial'
                        call mpi_abort (mpi_comm_world, 1, error)
                     endif
                  endif

                  j_block_n_p = j_block_n
                  if (any (bnd_exterior(:, :, 2) .eq. 0) .and. (j_block_n .eq. 0 .or. j_block_n .eq. ini_block_y + 1)) then
                     if (j_block_n .eq. 0) then
                        j_block_n_p = ini_block_y
                     elseif (j_block_n .eq. ini_block_y + 1) then
                        j_block_n_p = 1
                     else
                        write (*, *) 'Periodic boundary error: model_initial'
                        call mpi_abort (mpi_comm_world, 1, error)
                     endif
                  endif

                  k_block_n_p = k_block_n
                  if (any (bnd_exterior(:, :, 3) .eq. 0) .and. (k_block_n .eq. 0 .or. k_block_n .eq. ini_block_z + 1)) then
                     if (k_block_n .eq. 0) then
                        k_block_n_p = ini_block_z
                     elseif (k_block_n .eq. ini_block_z + 1) then
                        k_block_n_p = 1
                     else
                        write (*, *) 'Periodic boundary error: model_initial'
                        call mpi_abort (mpi_comm_world, 1, error)
                     endif
                  endif
                  if (i_block_n_p .ge. 1 .and. i_block_n_p .le. ini_block_x .and.                                                  &
                      j_block_n_p .ge. 1 .and. j_block_n_p .le. ini_block_y .and.                                                  &
                      k_block_n_p .ge. 1 .and. k_block_n_p .le. ini_block_z) then
                     child%neigh(i_neigh, j_neigh, k_neigh)%pblock => base%child(i_block_n_p, j_block_n_p, k_block_n_p)%pblock
                  endif
               endif
            enddo
         enddo
      enddo

   enddo
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  distribution blocks for different cpu
!----------------------------------------------------------------------------------------------------------------------------------|
   blocks_front = 0
   do i_cpu = 1, cpu_id
      blocks_front = blocks_front + n_blocks(i_cpu)
   enddo

   tmp_start = blocks_front + 1
   tmp_end = blocks_front + n_blocks(cpu_id + 1)

   call mpi_allgather(tmp_start, 1, mpi_integer4, hilbert_start, 1, mpi_integer4, mpi_comm_world, error)
   call mpi_allgather(tmp_end, 1, mpi_integer4, hilbert_end, 1, mpi_integer4, mpi_comm_world, error)

   current => base%next
   do while (current%lv .ne. 0)
      if (current%hilbert_number .eq. hilbert_start(cpu_id_p1)) then
         current%last => base
         base%next => current
      endif
      if (current%hilbert_number .eq. hilbert_end(cpu_id_p1)) then
         base%last => current
         current%next => base
      endif
      current => current%next
   enddo

   if (cpu_number .gt. n_blocks_all) then
      if (cpu_id_p1 .gt. n_blocks_all) then
         base%next => base
         base%last => base
      endif
   endif

   current => base%framework_next
   do while (current%lv .ne. 0)
      if (current%hilbert_number .lt. hilbert_start(cpu_id_p1) .or. current%hilbert_number .gt. hilbert_end(cpu_id_p1)) then
         current%last => null()
         current%next => null()
      endif
      current => current%framework_next
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  allocatememory and calculate the maximum speed
!----------------------------------------------------------------------------------------------------------------------------------|
   max_speed = -1.0d10
   current => base%framework_next
   do while (current%lv .ne. 0)

      if (current%cpu_id .eq. cpu_id) then
         allocate(current%ro(nx_block, ny_block, nz_block), current%mx(nx_block, ny_block, nz_block),                              &
                  current%my(nx_block, ny_block, nz_block), current%mz(nx_block, ny_block, nz_block),                              &
                  current%bx(nx_block, ny_block, nz_block), current%by(nx_block, ny_block, nz_block),                              &
                  current%bz(nx_block, ny_block, nz_block), current%en(nx_block, ny_block, nz_block),                              &
                  current%gx(nx_block, ny_block, nz_block), current%gy(nx_block, ny_block, nz_block),                              &
                  current%gz(nx_block, ny_block, nz_block), current%po(nx_block, ny_block, nz_block))
         allocate(current%fxro_bnd(ny_block, nz_block, n_neighs), current%fxmx_bnd(ny_block, nz_block, n_neighs),                  &
                  current%fxmy_bnd(ny_block, nz_block, n_neighs), current%fxmz_bnd(ny_block, nz_block, n_neighs),                  &
                  current%fxbx_bnd(ny_block, nz_block, n_neighs), current%fxby_bnd(ny_block, nz_block, n_neighs),                  &
                  current%fxbz_bnd(ny_block, nz_block, n_neighs), current%fxen_bnd(ny_block, nz_block, n_neighs),                  &
                  current%fxpo_bnd(ny_block, nz_block, n_neighs),                                                                  &
                  current%fyro_bnd(nx_block, nz_block, n_neighs), current%fymx_bnd(nx_block, nz_block, n_neighs),                  &
                  current%fymy_bnd(nx_block, nz_block, n_neighs), current%fymz_bnd(nx_block, nz_block, n_neighs),                  &
                  current%fybx_bnd(nx_block, nz_block, n_neighs), current%fyby_bnd(nx_block, nz_block, n_neighs),                  &
                  current%fybz_bnd(nx_block, nz_block, n_neighs), current%fyen_bnd(nx_block, nz_block, n_neighs),                  &
                  current%fypo_bnd(nx_block, nz_block, n_neighs),                                                                  &
                  current%fzro_bnd(nx_block, ny_block, n_neighs), current%fzmx_bnd(nx_block, ny_block, n_neighs),                  &
                  current%fzmy_bnd(nx_block, ny_block, n_neighs), current%fzmz_bnd(nx_block, ny_block, n_neighs),                  &
                  current%fzbx_bnd(nx_block, ny_block, n_neighs), current%fzby_bnd(nx_block, ny_block, n_neighs),                  &
                  current%fzbz_bnd(nx_block, ny_block, n_neighs), current%fzen_bnd(nx_block, ny_block, n_neighs),                  &
                  current%fzpo_bnd(nx_block, ny_block, n_neighs))

         call model_variables(current%ro, current%mx, current%my, current%mz, current%bx, current%by, current%bz, current%en,      &
            current%gx, current%gy, current%gz, current%po, current%x, current%y, current%z, nx_block, ny_block, nz_block)

         if (ini_value_save_flag .ne. 0) then
            allocate(current%ro0(nx_block, ny_block, nz_block), current%mx0(nx_block, ny_block, ny_block),                         &
                     current%my0(nx_block, ny_block, ny_block), current%mz0(nx_block, ny_block, ny_block),                         &
                     current%bx0(nx_block, ny_block, ny_block), current%by0(nx_block, ny_block, ny_block),                         &
                     current%bz0(nx_block, ny_block, ny_block), current%en0(nx_block, ny_block, ny_block))

            do k = 1, nz_block
            do j = 1, ny_block
            do i = 1, nx_block
               current%ro0(i, j, k) = current%ro(i, j, k)
               current%mx0(i, j, k) = current%mx(i, j, k)
               current%my0(i, j, k) = current%my(i, j, k)
               current%mz0(i, j, k) = current%mz(i, j, k)
               current%bx0(i, j, k) = current%bx(i, j, k)
               current%by0(i, j, k) = current%by(i, j, k)
               current%bz0(i, j, k) = current%bz(i, j, k)
               current%en0(i, j, k) = current%en(i, j, k)
            enddo
            enddo
            enddo
         endif

         current%var_allocate_flag = .true.

         do k = 1, nz_block
         do j = 1, ny_block
         do i = 1, nx_block
            vx = current%mx(i, j, k) / current%ro(i, j, k)
            vy = current%my(i, j, k) / current%ro(i, j, k)
            vz = current%mz(i, j, k) / current%ro(i, j, k)
            b2 = current%bx(i, j, k) * current%bx(i, j, k) +                                                                       &
                 current%by(i, j, k) * current%by(i, j, k) +                                                                       &
                 current%bz(i, j, k) * current%bz(i, j, k)
            v2 = vx * vx + vy * vy + vz * vz
            pr = (current%en(i, j, k) - v2 * current%ro(i, j, k) / 2.0d0 - b2 / 2.0d0) * gmm1
            max_te_ro = max(max_te_ro, (pr / current%ro(i, j, k)) ** (5.0d0 / 2.0d0) / current%ro(i, j, k))
            c2 = gm * pr / current%ro(i, j, k)
            s2 = c2 + b2 / current%ro(i, j, k)
            cax = abs(current%bx(i, j, k)) / sqrt(current%ro(i, j, k))
            cay = abs(current%by(i, j, k)) / sqrt(current%ro(i, j, k))
            caz = abs(current%bz(i, j, k)) / sqrt(current%ro(i, j, k))
            cfx = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * cax * cax)) / 2.0d0)
            cfy = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * cay * cay)) / 2.0d0)
            cfz = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * caz * caz)) / 2.0d0)
            max_speed = maxval((/ max_speed, cfx + abs(vx), cfy + abs(vy), cfz + abs(vz) /))
         enddo
         enddo
         enddo
      endif

      current => current%framework_next
   enddo

!   current => base%framework_next
!   do while (current%lv .ne. 0)
!      if (cpu_id .eq. 0) print*, current%hilbert_number, current%cpu_id, cpu_id, allocated (current%ro)
!      current => current%framework_next
!   enddo
!pause
!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine model_initial
